CC written by M. Widom CMU march 2004
c     compile with g77 -llapack
c     input is:(1) TM (2) List of T values
c     TM format:dimension followed by list of elements 1/line
c     element format: # of terms, coeff and exponent for each term.
c     output is: T, internal energy and heat capacity (/kB)
      implicit none
      real*8 M(217,217),T,ev,u,c,kB,s
      integer coeff(197),nk(217,217),nt,dim,ipr
      real*8 expon(197),dT,T1,T2,T3,ev1,ev2,ev3,xnorm
      common/mat/nk,coeff,dim
      common/expn/expon
      dT=1.E-3
      read(*,*)xnorm
      call ReadMat
      open(9,file='tmevec.dat')
      do nt=1,10000
         read(*,*,end=1) T
         T1=T*(1.-dT)
         T2=T
         T3=T*(1.+dT)
         call getEV(T1,ev1,0)
         call getEV(T2,ev2,1)
         call getEV(T3,ev3,0)
         call MakeThermo(T1,T2,T3,ev1,ev2,ev3,u,c,s)
         write(*,900) T,c/xnorm,u/xnorm,T*s/xnorm
      enddo
 1    continue
 900  format(f10.2,3f15.6)
      stop
      end

      subroutine ReadMat
      integer dim,nk(217,217),i,j,k,coeff(217,217,197)
      real*8 expon(217,217,197)
      common/mat/nk,coeff,dim
      common/expn/expon
      read(*,*) dim
      if(dim.gt.217) then
         write(*,*) 'Index out of bounds',dim
         stop
      endif
      do i=1,dim
         do j=1,dim
            read(*,*) nk(i,j),(coeff(i,j,k),expon(i,j,k),k=1,nk(i,j))
            if(nk(i,j).lt.0.or.nk(i,j).gt.197) then
               write(*,*) 'Index out of bounds',i,j,nk(i,j)
               stop
            endif
         enddo
      enddo
      return
      end

      subroutine getEV(T,ev,ipr)
      implicit none
      real*8 T,ev
      real*8 kB
      integer dim,nk(217,217),ipr,i,j,k,coeff(217,217,197)
      real*8 expon(217,217,197)
      common/mat/nk,coeff,dim
      common/expn/expon
      integer LDVL,LDVR,LWRK,INFO,ii,io
      real*8 WRK(4000),VR(217,217),VL(217,217),WR(217),WI(217)
      real*8 M(217,217),elm
      character*1 job,uplo
      kB=1./11604.
      LWRK=4000
      LDVL=1
      LDVR=1
      do i=1,dim
         do j=1,dim
            elm=0.0
            do k=1,nk(i,j)
               elm=elm+coeff(i,j,k)*dexp(-expon(i,j,k)/(kB*T))
            enddo
            M(i,j)=elm
         enddo
      enddo
      job='N'
      uplo='N'
      call DGEEV(job,uplo,dim,M,217,WR,WI,VL,LDVL,VR,LDVR,WRK,LWRK,INFO)
c      if(ipr.eq.1) then
c         ev=-100000.0
c         do i=1,dim
c            if(WR(i).gt.ev) then
c               ev=WR(i)
c               ii=i
c            endif
c         enddo
c         write(*,2791)T,ii,ev
c         do i=1,dim
c            write(9,8433)T,(vr(ii,j),j=1,dim)
c         enddo
c 2791    format(f10.1,i7,f12.6)
c 8433    format(f8.1,1000f10.5)
c      endif
      return
      end

      subroutine MakeThermo(T1,T2,T3,ev1,ev2,ev3,u,c,s)
      real*8 T1,T2,T3,ev1,ev2,ev3,u,c,s
      real*8 f1,f2,f3,u15,u25,kB,T15,T25
      kB=1./11604.
      T15=(T1+T2)/2.
      T25=(T2+T3)/2.
      f1=-kB*T1*dlog(ev1)
      f2=-kB*T2*dlog(ev2)
      f3=-kB*T3*dlog(ev3)
      u15=-T15*T15*(f2/T2-f1/T1)/(T2-T1)
      u25=-T25*T25*(f3/T3-f2/T2)/(T3-T2)
      u=(u15+u25)/2.
      c=(u25-u15)/(T25-T15)/kB
      s=(u-f2)/T2
      return
      end
